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Abstract 

The properties of UO2 result from rich /-electron physics, including electronic Coulomb inter- 
actions, spin-orbit and crystal field effects, as well as inter-ionic multipolar coupling. We present 
a comprehensive theoretical study of the electronic structure of UO2 using a combined application 
of self-consistent DFT+£7 calculations and a model Hamiltonian. The T§ ground state of U 4+ and 
the energies of crystal field excitations T§ — >► I^^i are reproduced in very good agreement with 
experiment. We also investigate competing non-collinear magnetic structures and confirm 3-k as 
the T = 0K ground state magnetic structure of UO2. 

PACS numbers: 71.27.+a, 71.15.Mb, 71.70.Ch 



I. INTRODUCTION 



Uranium dioxide is an important and interesting material from both technological and 
scientific perspectives. During the past half-century, the electronic structure of UO2 has 
been thoroughly characterized by various experiments^! (f or a recent review, see Ref. [T3l) . 
UO2 is a semiconductor with a 2 eV band gap 1 and localized 5/ 2 -electrons that retain strong 
atomic-like properties. Due to significant Coulomb interactions and spin-orbit (SO) effects, 
the ground state of a free U 4+ ion is the 3 i^4 nonet [see Fig. [jja)]. When the crystal-field 
(CF) of UCVs fluorite structure is considered, 3 #4 is split into the ground state T 5 triplet and 
the excited T 3 doublet, T 4 triplet, and Ti singlet, all approximately 0.15 eV above [see 
Fig. [1(b)]. When cooled below T N = 30.8 K, U0 2 under goes a first-order phase transition 
from a paramagnetic to a transverse type-I antiferromagnetic (AF) phase 2 , which exhibits 
a Jahn- Teller (JT) distortion of the oxygen cage 5 . Experimental studies now converge on 
the view that the non-collinear magnetic structure and the oxygen distortion are of the 3-k 
typ^^U^, i.e., the moment and lattice distortion are both along the (111) direction [see 
Fig.Qc)], instead of the previously proposed 1-k ((001))Pand 2-k ((110) jl 3 structures. 




free U 4+ (5f 2 )ion cubic crystal field coupling with whole lattice 

FIG. 1. Schematics of the 5/ 2 ground states and level splitting, in decreasing interaction strength, 
of (a) free U 4+ ion, (b) cubic CF, and (c) ordered 3-k non-collinear magnetic structure of bulk UO2: 
left, direction of magnetic moments on uranium, designated by large arrows; right, distortion of 
oxygen (small arrows) around a central U atom. 

On the theory side, the CF model of Rahman and Runcimaif^ correctly predicted the T 5 
ground state of UO2 (Fig. [2]). Recent CF calculations have obtained quantitative agreement 
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with experimental excitation spectra by fitting model parameters to the measured dataF^, 
by adding corrections to the point charge model 17 , or by extrapolating from the fitted values 
for other actinide dioxide^. Models of magnetism in UO2, pioneered by the work of Alleif^, 
have explored the delicate interplay between multipolar and Jahn- Teller effect^. 

First-principles calculations have to go beyond the local-density or generalized-gradient 
approximations (LDA/GGA) to the density-functional theory (DFT) to correctly reproduce 
the insulating character of UO2. Existence of an energy gap was demonstrated in Ref. [T9l 
and |20] using the hybrid functional methocP^, in Ref. [22] using the self-interaction-corrected 
LDA 23 , and in Ref. El using the DFT+U method^. CF splitting in actinide compounds has 
been computed by using constrained / states without full self-consistency 26 or by analyzing 
band positions obtained from LDA/GGA calculation^. The 3-k structure of UO2 was 
studied by Laskowski et al using DFT+[/ 28 , but their results showed anomalous dependence 
on the U parameter and both the calculated oxygen distortions and energy differences were 
about an order of magnitude too large. Furthermore, the 3-k state was only stable with 
large U values and a formulation of DFT+U that is usually only applied to metals. A 
first-principles framework for self-consistently and accurately accounting for all the different 
energy scales in Fig. [I] does not yet exist. 

In this paper, we present a unified DFT-based framework for calculating the electronic 
spectra, magnetism and lattice distortions in UO2. Explicit /-/ interactions and CF ef- 
fects are treated using a model Hamiltonian with parameters derived from self-consistent 
DFT+U calculations. The ground state wave functions that are obtained by diagonalizing 
this Hamiltonian are used to set up initial conditions for self-consistent DFT+U calculations 
of magnetism and lattice relaxations. Our approach allows us to accurately reproduce all 
the different energy scales in Fig. [TJ including the T 5 ground state, F 3 ^i excited states, as 
well as the energetics of competing magnetic structures, including 3-k, and their associated 
lattice distortions, all within a unified self-consistent framework. 

Before moving on to the details, we would like to stress that extra care should be taken 
in first-principles calculations of / electrons. Several challenges are encountered in DFT 
calculations of UO2 (and other actinide compounds in general). First, strong /-/ inter- 
actions and a weak crystal field result in an inherently complicated many-body problem. 
For instance, since 5/ 2 electrons hybridize weakly with the O 2p bands and remain well 
localized, their true wavefunctions are in general multi-determinantal (see below for further 
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discussions). Second, the higly localized nature of / electrons tends to magnify the inaccu- 
racies of approximate exchange-correlation functionals. We have previously shown that the 
self-interaction (SI) error of /-electrons is highly sensitive to the occupied orbital, and its 
removal is non-trivial in both the DFT+U and hybrid-functional method^. Therefore, an 
improved version of DFT+c^O is required to remove such errors (~ 0.1 eV) and access weak 
CF effects. Third, the existence of a multitude of f-states often leaves DFT+U calculations 
trapped in local minima, leading to difficulties in reproducibly finding the correct electronic 
ground state. As a consequence, it is not uncommon for different authors to find inconsistent 
and hard-to-interpret results with large errors (~ 1 eV or even larger) even when using iden- 
tical electronic-structure methods (see Refs. I29H361 and references therein). It is likely that 
this issue contributed to the failure of previous studies to reliably examine the 3-k structure 
of UO2. 28 Previously, we have shown that the local minima issue is also present in hybrid 
functional calculations 29 . In this paper, we show that the multiple minima, corresponding 
to different orbital states, contain valuable physical information about /-electrons that can 
be used to help find the true ground state and excitation spectra. 
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FIG. 2. Low energy f n eigenstates of Hamiltonian [2] a) n = 1: Tg ground states and doublet of 
the j = 5/2 sextet b) n — 2: T§ ground states and the excited ^z^i of 3 i^4 



II. METHOD 

A. LDA+J7 calculations 

All DFT calculations were carried out using the VASP codeP, GGA-PAW potentiakP, a 
cutoff energy of 450 eV, and without any symmetry constraints to ailow symmetry-breaking 
solutions. Crystal field caiculations were performed in the primitive cell of one UO2 formula 
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unit with a 6 x 6 x 6 fc-point grid. The lattice and ionic positions were frozen at the 
experimental fluorite structure for crystal field calculations. These calculations, as discussed 
in sectioi jTTl A[ are ferromagnetic with one uranium ion per cell. Magnetic structures were 



calculated in the fee supercell (4 formula units) using a 4 x 4 x 4 grid, first without and 
then with full relaxation. Spin-orbit coupling was self-consistently incorporated for realistic 
comparison with experiment. 

To remove the orbital-dependent components of self-interaction errors (SIE) of f- 
electrons, we use a formulation of the LDA+U method^ by modifying only the exchange 
term, rather than both Hartree and exchange, of the LDA: 



where the orbital-dependent Hartree- Fock exchange Ex contains a term that approximately 
cancels the on-site SIE in the Hartree energy of localized /-electrons; the remainder of the 
LDA Hartree energy is exact by definition and therefore left unmodified in our approach. 
The exchange double-counting term E^cx accounts for the LDA exchange energy and is given 
by a linear combination controlled by the c parameter of the exchange double-counting in 
the LiechtensteiiP^ scheme and the on-site local-spin-density (LSD) exchange, conceptually 
similar to hybrid functional approaches and serves the purpose of subtracting the orbital- 
dependence of the LDA exchange energy. As a result, eq. [T| is self-interaction free to high 
accuracy. 

There is only one adjustable parameter, U ) in our approach, and the other parameters 
J and c can be determined at given U. As done in Ref. [29l we choose up to seven f 2 SSD 
states of the U 4+ ion that are analytically degenerate without considering spin-orbit, and 
calculate these states' total energy dependence on J and c. As shown in Fig. |3j optimal 
values of J=0.6 eV and c = 0.5 are obtained at U =6 eV that minimize the energy difference, 
i.e. the orbital-dependent self-interaction error. These J and c values are used throughout 
the paper. We use U =6 eV in this paper and discuss the dependence of the results on U in 



section III B In the rest of the paper spin-orbit is included. 
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FIG. 3. LDA+U energy of U 4+ ion as function of J in different orbitals. SOC is not included. 
B. On-site model Hamiltonian for / 

We consider the following single-ion model for /-electrons: 

n 

H = J2(fi + (ii-Si) + V ee , (2) 

1=1 

where the summation runs over n electrons for the one-body terms of cubic CF, /, and SO 
coupling of strength (. The electronic interaction V^e is parametrized by Slater's integrals F k 
{k = 0, 2, 4, 6)3^. The matrix elements of / between the basis states indexed by projections 
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of orbital (m) and spin (a) momenta are given by 
(ma\f\mW)=6 fTfT , I Y> ^ 



V 4 (Y 4 
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32 x /-V 6 (Y m - V145RF 64 ) 



Yi,dn, 



(3) 



where V^q are cubic CF parameters^ and are complex spherical harmonics. To study 
the magnetic properties, an infinitesimal magnetization field B (B 0) is applied: 

n 

H ' = h ~Y, b ' (aJi + gsSi)fi B /fr, (4) 
i=i 

where ql = 1 and gs ~ 2 are the orbital and spin g-f actors, respectively. 
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TABLE I. For different magnetization directions B, the ground states quartet (f 1 ) and T$ 
triplet (/ 2 ), their spin and total magnetic moment in and, for / 2 , the dominant determinants 
in the corresponding f 1 basis. 

We first discuss the general properties of solutions to Eqs. (J2])- (J4]) using the model param- 



eters derived from DFT+£7 calculations (which will be discussed in detail in section III A). 
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For n = 1, 14 eigenstates are obtained, the lowest being the T 8 quartet (Fig. [2] and Table [T]). 
For n > 1, the Hamiltonian in Eq. ^ can be diagonalized via configuration interaction of 
C* 4 f n single Slater determinants (SSD) based on the f 1 eigenstates. The T 5 ground states 
of f 2 are shown in Table [TJ together with their dominant determinants, designated as 
using the indices of f 1 states in the left column of Table |lj For magnetic moment along each 
of the [001], [110], and [111] directions, the T 8 quartet of f 1 includes states 1,4 (2,3) with 
larger (smaller) spin and orbital magnetic momenta, while the T 5 triplet of f 2 consists of 
states a and c with = 2.06 \ib and one dominant determinant [(1, 2) or (3, 4)], as well as 
a non-magnetic state b dominated by two determinants (right column of Table [T]). Note that 
the observed moment of the ordered state is smaller at 1.75 /i^F. The moment fi = 2.06/i^ of 
the T 5 (a, c) states is slightly larger than the saturated 2fi B characteristic of the 3 i?4 multiplet 
because exchange and SOC interactions are all of comparable strength and other multiplets 
slightly mix into the ground state and increase the effective momentPH In general, all the 
f 2 eigenstates, including T^^i (Fig. [2]) are composed of multiple determinants. 



C. Model parameters from LDA+£7 

The parameters for the model Hamiltonian in Eq. ^ are obtained by analyzing the 
total energies and / wavefunctions calculated with LDA+J7. In this procedure, many self- 
consistent LDA+U calculations are first carried out, yielding solutions that are in general 
local minima rather than the global minimum of UO2. Next, we extract the f 2 Kohn-Sham 
wavefunction |^/), SSD by construction, from each solution, and compute the expected 
energy according to Eq. ^ 

(^ f \H \^ f ) = x,V 4 + x 2 V 6 + x 3 ( + X4F 2 + U, (5) 

where a^'s represent the solution-dependent coefficient associated with model parame- 
ters. Since the F k {k = 2,4,6) contributions of V ee are heavily correlated^!, the following 
approximatioiP^ has been adopted in Eq. ^ : 

F 2 = F 4 /0.668 = F 6 /0.494, (6) 

eliminating model parameters F 4 and F 6 . Finally, expectation values of H (Eq. [5]) of 
the obtained solutions are fitted to the corresponding DFT+U total energies, yielding self- 
consistent ab initio values of F 2 , ( and the crystal field parameters V4 and Vq. We use the 
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simple least-square method to perform the linear fitting. Here U in eq. ^ can be regarded 
as a constant in the fitting and bears no direct physical meaning. 



III. RESULTS AND DISCUSSIONS 



A. Crystal field ground states and excitations 



We carried out a series of 50 different self-consistent calculations with randomly initialized 
f 2 states. Due to the existence of multiple local minima in DFT+£7, these calculations 
resulted in a range of energies spread over almost 2 eV [filled circles in Fig. |2|(a)]. It is seen 
that random wavefunction initialization has generated only one low-energy solution, while 
the remaining runs were trapped in met ast able high-energy states. 
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FIG. 4. a) Fitting of BFT+U 29 energy to Eq.[2]for 50 runs with random initial states (filled circles) 
and 15 states with initial states constructed from f 1 solutions (open circles), b) Predicted f 2 CF 
levels T5 5 3 5 4 5 i compared to measured CF splitting. 



Obtained from the fitting procedure outlined in section II C , model parameters are ap- 
plied in Eq. ^ to construct f 1 eigenstates and subsequently determine f 2 states by direct 
diagonalization within the subspace of SSDs formed from f 1 eigenstates. To further improve 
the quality of our fit and provide data points in the low-energy region that was poorly rep- 
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resented in the randomly initialized sample [filled circles in Fig. [4](a)], we self-consistently 
calculate the DFT+U energies of additional 15 two-electron SSD states that involve the 6 
low energy f 1 orbitals (T 8 and F 7 ) for B//[001]; these points are shown as open circles in 
Fig. [4^. The (1,2) and (3,4) states, which dominate the T 5a and r 5c ground states, are 
also found to have the lowest energies in self-consistent DFT+U calculations, demonstrat- 
ing that our method can reliably locate the electronic ground state. The other ground state 
in Table |TJ T 5 5, has two dominant determinants and is not directly accessible in DFT+U. 
Therefore, the data flow between the model Hamiltonian and DFT+U calculations is bi- 
directional: DFT+U provides model parameters, while the model guides the DFT+U to 
the ground state and gives access to multi-determinant states. 
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The final fitted parameters are shown in Table [TTJ Compared to the values obtained by 
fitting the spectra of free ion^H, the ionic parameters F k and ( in the UO2 solid are somewhat 
suppressed due to hybridization and screening effects. The calculated cubic CF parameters 
V4 and Vq are slightly smaller than those fitted to experimental data or extrapolated from 
other actinide oxided^^. The Hamiltonian in Eq. ^ can now be diagonalized. The 
predicted energies of the three lowest excited CF levels F 3 ^i are in reasonable agreement 
with experimentP with errors of approximately 10-20 meV [Fig. [l](b)]. A notable deviation 
is overestimation of the splitting between these levels. 

Finally, we note that the input J = 0.6 eV used in our DFT+U calculation differs from 



the fitted value of J' = (286F 2 + 195F 4 + 250F 6 )/6435 in Table (n) This is because the 
role of the former is to minimize the SI error in DFT+£7, while the latter represents on-site 
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exchange, and some difference between them is expected when used with an approximate xc 
functional. A perfect xc functional would make the input U or J unnecessary and predict 
physically meaningful output J f or F k . 



B. Dependence on input U 

To illustrate the effect of the only adjustable variable in our approach, U : the same 
calculations were repeated using U = 4.5 eV. As shown in Table and Fig.[I](b), the results 
change only slightly and remain in good agreement with experiment. Note that when J7, 
which controls the degree of electron localization, is decreased, the ionic parameters F k 
and ( also decrease, i.e. away from the free ion values, while the CF parameters increase, 
suggesting that the /-electrons become more delocalized. Such a picture of opposite influence 
of electron localization on free ion and CF parameters is consistent with the observed trend 
that increase of the CF interaction results in a reduction in the free-ion parameters for the 
same ion in different chemical environment!^. 



C. Magnetic properties 

Finally, we discuss the effects of inter-ionic interactions and magnetic properties of UO2. 
Various magnetic structures within a cubic supercell of four formula units were calculated, 
first without and then with lattice relaxation. The previous approximation of representing 



the T 5 ground state with the (1,2), (3,4) SSDs was adopted. Table III shows the energies 
E s (static lattice), E r (after full relaxation), and the total magnetic moments \i assuming 
ferromagnetic (FM) and type-A anti-ferromagnetic (AAF) configurations along the [001] 
(reference), [110], and [111] directions. In all cases, the calculated \i of ~ 2.1-2.2 \i B is 
close to the saturated value 2.06 \ib in Tableland larger than the measured 1.74 /x^ 2 . The 
reduction of the ordered moment is a topic of considerable interest, but the mechanisms, 
such as the dynamical Jahn- Teller effect 46 , are beyond the scope of this work. Given that 
the (l,2)/(3,4) states are the CF ground states within our computational approach and 
the calculated moment is not too far from 1.74 we continue to use these settings for 
non-collinear magnetic calculations. 



Table [TTT] shows that for each magnetization direction, the AAF configuration is always 
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lower than the FM configuration, in agreement with experiment. The energy differences 
E S (FM) — E S (AAF) are in the 4-6 meV range, suggesting that multipolar interactions and 
anisotropy are weak. These weak, mostly isotropic interactions underlie the success of the 



simple fitting procedure of Section II C (root-mean-square error = 35 meV) in FM config- 
urations of 15 (001) solutions and 50 solutions with random magnetic moment directions. 
We also find that the different magnetization directions differ in energy by less than 9 meV, 
suggesting that our procedure for removing the orbital-depedent SI error is highly accurate. 
Indeed, the f 2 wavefunctions differ considerably for the three principle directions (see Table 
[]]), which would result in SI errors of 0.1 ^ 0.2 eV using the unmodified DFT+U approach. 
Nevertheless, we take additional care to remove any remaining SI errors, however small they 
they appear to be, by subtracting a reference energy for each principle magnetization 
direction d = (001), (110), or (111): E d = [E a (FM) + E s (AAF)}/2. With this correction, the 
AAF configurations in the three directions, as well as the 3-k structure (magnetic moment 
along (111), see Fig. [TJ^), are essentially degenerate before lattice relaxation. 

The magnetic transition temperature of antiferromagnetic UO2 is estimated with a clas- 
sical Heisenberg model on an fee lattice: 

U = -JH^Tsi ' Sj, (7) 

where the summation is over all nearest neighbor sites (ij) with unit spin s. The FM/AAF 
energy difference per UO2 is AE = 6J# — (— 2J#) = 8J#. As shown in Table [TTT| AE ~ 6 
meV, corresponding to T N = 3.18 Jh /ks^ 7 ) or about 28 K, in excellent agreement with the 
experimental value 30.8 K. 2 

After relaxation without symmetry constraints, we obtain the energies listed in the right 



column of Table III, The computed moments \i become slightly larger than the static values. 
The relaxation energy E r and the corrected E r — E^ are large (> 50 meV) due to differences 
between the static (fixed to experimental a = 5.47 A) and relaxed lattice parameters. The 
energy differences between competing magnetic configurations increases to ~ 10 meV; these 
values are consistent with the Neel temperature of T N = 30.8 K 2 . The relaxed structures 
with fi/ /(111) are clearly more stable than (001) and (110). We have enumerated all AF 
(111) structures within the fee unit cell and found that the 3-k structure indeed has the 
lowest energy. The associated oxygen lattice distortion (amplitude 0.024 A) is also of the 
3-k type, though slightly larger than the measured 0.014 A 5 . 
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TABLE III. Energy (in meV per UO2) for different magnetic structures, without and with ionic 
relaxation. 

IV. CONCLUSIONS 

In summary, we have studied the electronic structure of UO2 using an aspherical-self- 
interaction free DFT+£7 method coupled with a model Hamiltonian. The T 5 crystal field 
ground states, as well as the CF excitation energies are reproduced in good agreement with 
experiment. Various magnetic structures are investigated with careful initialization of the 
orbital and magnetic states. The inter-ionic interactions are found to be weak and largely 
isotropic. When self-interaction errors are accounted for, the 3-k structure is essentially 
degenerate with other antiferromagnetic configurations and becomes the ground state only 
when lattice relaxations are considered. Our work demonstrates the usefulness of electronic 
structure calculations for /-compounds with proper treatment of self-interaction errors and 
multiple self-consistent local minima corresponding to different orbital states; this approach 
can be readily applied to defect supercells and other /-compounds. 
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